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Robust edge states and non-Abelian excitations are the trademark of topological states of matter, 
with promising applications such as "topologically protected" quantum memory and computing. 
While so far topological phases have been exclusively discussed in a Hamiltonian context, we show 
that such phases and the associated topological protection and phenomena also emerge in open 
quantum systems with engineered dissipation. The specific system studied here is a quantum wire 
of spinless atomic fermions in an optical lattice coupled to a bath. The key feature of the dissipa- 
tive dynamics described by a Lindblad master equation is the existence of Majorana edge modes, 
representing a non-local decoherence free subspace. The isolation of the edge states is enforced by 
a dissipative gap in the p-wave paired bulk of the wire. We describe dissipative non-Abelian braid- 
ing operations within the Majorana subspace, and we illustrate the insensitivity to imperfections. 
Topological protection is granted by a nontrivial winding number of the system density matrix. 



Topological properties can protect quantum systems 
from microscopic details and imperfections. In condensed 
matter physics this is illustrated by the seminal exam- 
ples of the quantum Hall effect and the recently discov- 
ered topological insulators [iHSj . The ground state of the 
Hamiltonian of such systems is characterized by nonzero 
values of topological invariants which imply the existence 
of robust edge states in interfaces to topologically trivial 
phases. Due to their topological origin, these modes are 
immune against a wide class of perturbations. 

The conceptually simplest example illustrating these 
phenomena is Kitaev's quantum wire representing a 
topological superconducting state supporting Majorana 
fermions as edge states [7]. The pair of Majorana edge 
modes represents a nonlocal fermion which is a promising 
building block to encode topological qubits [HHTO] . Simi- 
lar to the Majorana excitations near vortices of a Px+Wy 
superconductor [TTl [12] , they show nonabelian exchange 
statistics when braided in ID wire networks [10]. 

Remarkably, the above described topological features 
and phenomena not only occur as properties of Hamil- 
tonians, but appear also in driven dissipative quantum 
systems. Below we will develop such a topological pro- 
gram for a dissipative many-body system parallel to the 
Hamiltonian case. We will do this for a dissipative ver- 
sion of Kitaev's quantum wire. This represents the sim- 
plest instance exhibiting the key features such as dissipa- 
tion induced Majorana edge modes, decoupled from the 
dynamically created p-wave superfluid bulk by a dissi- 
pative gap. The dissipation induced topological order is 
generated in stationary states far away from thermody- 
namic equilibrium which are not necessarily pure, i.e. de- 
scribed in terms of a wave function only, and is reached 
exponentially fast from arbitrary initial states. This is 
in marked contrast to recent ideas of topological order 
in Hamiltonian systems under non-equilibrium periodic 
driving conditions [HI [T^ . Such a system can be realized 
with cold atoms in optical lattices, where the generation 



of topological order in the more conventional Hamilto- 
nian settings has been proposed recently in a variety of 
settings [T5HT9] . 

DISSIPATIVE EDGE MODES IN A FERMIONIC 
QUANTUM WIRE 

Our goal is to develop a master equation for a dissi- 
pative quantum wire which exhibits topological proper- 
ties including Majorana edge states. To illustrate the 
analogies and differences to the Hamiltonian case, and 
in particular to motivate our construction of the mas- 
ter equation with the topological states as dark steady 
states, we start by briefly summarizing Kitaev's model 
of the topological superconductor. 

Topological quantum wire - Kitaev considers spinless 
fermions a^, a] on a finite chain of N sites i described by 
a Hamiltonian 

N 

H = ^- Jaja^+i + (Aa^ai+i + h.c.) - /lajai , 

with a hopping term with amplitude J, a pairing term 
with order parameter A, and a chemical potential /i. The 
topologically non-trivial phase of the model is best illus- 
trated for the choice of parameters J = |A| and /i = 0, 
where the Hamiltonian simplifies to 

N-l N-1 

H = 2iJ ^ C2i C2i+i = 2 J ^ alcLi. (1) 

i=l i=l 

Here we have defined Majorana operators q as the 
quadrature components of the complex fermion operators 
o^i = I (ic2i-i +C2i) with properties = Cj,{cj,ci} = 
26 ji. The Hamiltonian is readily diagonalized in terms 
of fermionic Bogoliubov quasiparticle operators di = 
\ {c2i + ic2i+i), where importantly the pairing of Majo- 
ranas is from different physical sites. 
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FIG. 1. Schematic setup for the dissipative Majorana quan- 
tum wire, a) A reservoir represents a source and drain for 
the quantum wire which is coherent over each pair of sites. 
Independent of the initial condition, the bulk of the system 
is then cooled into a p-wave superfluid state by dissipatively 
establishing a pairing link for each two adjacent sites, b) Il- 
lustration of the stationary state, where in the Majorana basis 
of real fermions each physical site is split into two Majorana 
sites. In the bulk all Majorana modes from neighboring sites 
are paired (black links). For a finite wire two dissipative un- 
paired Majorana modes 7l, 7h appear at the edge as a highly 
nonlocal decoherence free subspace. They are isolated from 
the bulk by a dissipative gap. For a realization with cold 
atoms, see Fig. |2] 



The ground state satisfies the condition ai\G) = 
for all i. The bulk of the wire describes a fermionic 
p-wave superfluid with a bulk spectral gap^ which here 
equals the constant dispersion Ck = 2 J. For a finite 
wire, the absence of the term clj^cln = ic2ArCi indicates 
the existence of a two-dimensional zero energy non-local 
fermionic subspace spanned by \a) G {|0), |1) = ajy|0)}. 
While highly delocalized in terms of the original complex 
fermion, in the real Majorana basis the situation is de- 
scribed in terms of two Majorana edge modes = ci 
ilR = C2n) which are completely localized on the left- 
most (rightmost) Majorana site 1 (2A/'), describing "half 
a fermion each. These edge modes exist in the whole pa- 
rameter regime — 2J < /i < 2J, however leaking more 
and more strongly into the wire when approaching the 
critical values. Their existence is robust against pertur- 
bations such as disorder, which can be traced back to the 
bulk gap in connection with their topological origin [7 . 

Dissipative topological quantum wire - Consider the 
master equation for spinless fermions in a ID chain with 
N sites, 



dtp 



with p the system density operator. The two terms on 
the right hand side are a Hamiltonian and a dissipative 
term, respectively. Here we concentrate on purely dissi- 
pative dynamics which occurs at rate and set H = 0. 
We choose the Lindblad operators ji as the Bogoliubov 



operators defined above, 

ji = «i = = (3) 

These Lindblad operators are quasi-local superpositions 
of annihilation and creation operators (see Fig. [T]). Due 
to the fermionic nature of j^, this choice ensures that 
the bulk of the system "cools" under the above dynam- 
ics to the unique pure state ai\G) = 0, which by con- 
struction agrees with the p-wave superfluid ground state 
of the Hamiltonian ([T]). Following [20.-^22^ . this steady 
state is thus a many-body dark state of the Liouvillian, 
C{\G){G\) = 0. 

The approach to this steady state is governed by the 
damping spectrum of the Liouvillian C. Diagonality of 
C in the implies a flat damping spectrum hi = hi in 
analogy to the excitation spectrum of the Hamiltonian 
above. While the damping spectrum /^/c > is always 
positive semi-definite for fermions, this "dissipative gap" 
hzo = niin(/^/e) = implies exponentially fast approach of 
all observables to their steady state values. 

For a finite wire we find dissipative zero modes re- 
lated to the absence of the Lindblad operator a at- More 
precisely, there exists a subspace spanned by the edge- 
localized Majorana modes cln = ^ (iTl +7/?), with the 
above Fock basis \a) G {|0), |1)}, which is decoupled from 
dissipation, i.e. dtpa(3{t) = with pa(3 = {a\p\P). 

Edge modes as nonlocal decoherence free subspace - 
These dissipative edge modes are readily revealed in so- 
lutions of the master equation. Eq. Q is quadratic in 
the fermion operators, which implies solutions in terms 
of Gaussian density operators p{t) ~ exp [—^c^G{t)c]. 
Here we have defined a column vector c of the 2N Ma- 
jorana operators, and G is a real antisymmetric matrix 
related to the correlation matrix Tab{t) = ^([^0,^5]) = 
i[tanh(iG/2)]a6, which equally is real antisymmetric. 
Writing the Lindblad operators in the Majorana ba- 
sis, ji = ifc^jj — c^Vl^ such that the Liouvillian pa- 
rameters are encoded in a hermitian 2N x 2N matrix 
M = ^ • li (8) 4 5 ^his covariance matrix is seen to obey 
the dissipation-fluctuation equation [23j 



5tr = -{x,r}-F, 



(4) 



with real matrices X = 2ReM = and Y = 4ImM = 
—Y^. Physically, the matrix X describes a drift or 
damping, while the matrix Y is related to fluctuations 
in steady state, determined by {X, f } = —Y. Writing 
F = f + (5F, the approach to steady state is governed 
by dfST = — {X, ^F}, i.e., the eigenvalues of the positive 
semi-definite matrix X [24 give the damping spectrum. 
The "dark" nonlocal subspace of edge modes, decoupled 
from dissipation, is thus associated with the subspace of 
zero eigenvalues of the damping matrix X. In a spectral 
decomposition X = Ar|r)(r|, and identifying by greek 
subscripts the zero eigenvalues subspace, we can write by 
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FIG. 2. Microscopic implementation scheme for the Majo- 
rana LiouviUian Eqs. (2[3). The quantum wire is represented 
by the lower sites of an optical superlattice for spin polarized 
atomic fermions. They are coherently coupled to auxiliary 
upper sites by lasers with Rabi frequencies ±Q, alternating 
from site to site. Dissipation results from spontaneous Bogoli- 
ubov phonon emission via coupling of the system to a BEG 
reservoir (light grey). An edge can be created using single 
site addressability tools [261 127| , cutting off the lattice at an 
auxiliary site instead of a target system site. As shown in the 
implementation section, this setting reduces to the Lindblad 
operators ([3| at late times. 
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While the bulk {rs sector) damps out to the steady 
state by dissipative evolution, the density matrix in the 
edge mode subspace {aj3 sector) does not evolve, pre- 
serving its initial correlations. The coupling density ma- 
trix elements (mixed sectors) damp out according to 
^ri3{t) = e~^'^^rrf3{0); in the presence of a dissipative 
gap as in the example above, this fadeout of correlations 
is exponentially fast, leading to a dynamical decoupling 
of the edge subspace and the bulk. More generally, this 
structure of the master equation appears whenever there 
exists a basis in which each Lindblad operator is block 
diagonal with blocks associated to edge and bulk, and 
with vanishing entries in the edge block (see appendix). 

In summary, we arrive at the physical picture that 
dissipative evolution cools the bulk into a p-wave su- 
perfluid and thereby isolates the edge mode subspace, 
p{t oo) Pedge Pbulk, providing a highly nonlocal 
decoherence free subspace [25 . A physical implementa- 
tion of the master equation ([2| with cold atoms is out- 
lined schematically in Fig. [2] More details on the setup 
are given in the implementation section, following ideas 
of Ref. m. 



STABILITY OF EDGE MODE SUBSPACE 

Here we study the robustness of the edge mode sub- 
space against (i) global parameter changes in the Lind- 
blad operators, while preserving their translation invari- 
ance, and (ii) static disorder, which breaks this invari- 



ance. We consider two examples of quadratic master 
equations (|2| with Lindblad operators deviating from the 
ideal case Jsl), 



= ^ (sin 6>(ai - a^+i) +cos(9(a{ +a{^i)), (6) 
j^'' = 7^(sin6> (4 - tti+i) + cos 6 {ai + aJ^J), (7) 



where the ideal case corresponds to ^ = 7r/4. In the first 
case, the steady state of the bulk remains pure. In the 
second case, we find a mixed state while still preserving 
the properties of the edge subspace. As elaborated on 
in the appendix, this results from the fact that the first 
case (c) represents a canonical transformation up to nor- 
malization of the ideal Lindblad operators in momentum 
space, while the second one (n) is not (cf. Eq. (47) in 



the appendix). In this latter case, the steady state has no 
counterpart as a ground state of some Hamiltonian. This 
difference is illustrated in Figs. [2] a), b). There, we plot 
the purity spectrum spec(f^) in steady state (with the 
edge mode subspace initialized as pure); a pure state in 
the bulk is indicated by all eigenvalues of being equal 
to —1. Note that static disorder, implemented in terms 
of small random variations of the Lindblad parameters of 
range e <C ^ from site to site, makes the first case non- 
equivalent to a canonical transformation, and, therefore, 
degrades the purity of the steady state. 

While the purity spectrum is qualitatively different 
for both kinds of parameter deformations, the spectra 
of damping matrices, spec(X), are rather similar. The 
characteristic features are (i) the existence of a dissipative 
gap hzo = hzcosO^ closing at = 7r/2, and (ii) the exis- 
tence of two zero modes throughout the parameter space. 
The associated orthogonal eigenvectors vl,r describing 
the Majorana modes Jl^r = ^lr^ constructed 
explicitly (see appendix), with localization length given 
by hoc/a = (log Lln^lgos^ l)"' in both cases, in units of 
the lattice constant a. This shows the characteristic ex- 
ponential edge localization of Majorana modes close to 
the ideal case, while their extent becomes comparable to 
the system size close to the gap closing points (see Fig. 
|3|. Adding disorder modifies the bulk spectrum quanti- 
tatively of 0{e)^ while the zero mode subspace persists. 
In fact, the existence of two zero eigenvalues in the pres- 
ence of disorder can be established by explicit construc- 
tion (see appendix). 

In summary, a two-dimensional decoherence free sub- 
space exists for the entire parameter range notwithstand- 
ing disorder or the fact that the bulk steady state may 
be strongly mixed. A sensible notion of protection of the 
subspace in terms of dissipative isolation from the bulk, 
however, only exists sufficiently far away from the point 
where the damping gap closes. 
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FIG. 3. For a quantum wire with 50 lattice sites we plot 
the eigenvalue spectrum of the square steady state correlation 
matrix (Figs, a) and b)), and the eigenvalue spectrum of 
the damping matrix X (Figs, c) and d)) as a function of the 
eigenvalue index. In the last row, the amplitude moduli \vi\ 
for the left and right zero mode (edge mode) as a function of 
the lattice site index i are shown (Figs, e) and f)). The first 
and second column corresponds to quasi-canonical Q and 
non-canonical Lindblad operators ([7|, respectively. Results 
are shown for three angles: The ideal case — ^ Eq. ^ 
(green), — ^ (blue) and — ^ (red). Dissipative zero 
modes exist in all cases, but they become degenerate with 
the bulk modes in the third case (red) case where the bulk 
dissipative gap collapses. The closed circles show results for a 
globally fixed ^, while the open circles correspond to addition 
of local static disorder to the angles. 



ADIABATIC PARAMETER CHANGES AND 
DISSIPATIVE BRAIDING 

The above can be generalized to a dissipative quan- 
tum wire network of M finite chains, as counterpart of 
the Hamiltonian networks discussed by Alicea et al. [10] . 
This results in higher dimensional non-evolving sub- 
spaces of dimension 2M, again governed by the structure 
given by Eq. ([5|. As we will show below, this leads to the 
possibility of braiding of the dissipative Major ana edge 
modes by adiabatic parameter changes in the Liouvillian. 

We consider the time evolution of the density matrix 
in a CO- moving basis \a{t)) = U (t)|a(0)) which follows the 
decoherence free subspace of edge modes, i.e. preserves 
the property po^^ = 0. Demanding normalization of the 
instantaneous basis for all times, {b{t)\a{t)) = Sab^ this 
yields 



d 
dt 



P 



i[A,p]^^\a)pah{h\, 



(8) 



a,b 



with the unitary connection operator A = iU^U and 
pab = {a{t)\dtp\b{t)) the time evolution in the instanta- 



neous basis. The Heisenberg commutator clearly reflects 
the emergence of a gauge structure [28H32] in the density 
matrix formalism, which appears independently of what 
kind of dynamics - unitary or dissipative - generates the 
physical time evolution, represented by the second con- 
tribution to the above equation. The transformation ex- 
erted on the zero mode subspace of either Hamiltonian 
or Liouvillian with an initial condition pai3{0) is then 
given by pa/3{t) = with time-ordered 

V{t) = Texp (-i/o drA(r)) and A{t)^0 = 

Of central importance for such state transformations to 
work without losing the protected subspaces is adiabatic- 
ity of the parameter changes. Here, this is a requirement 
on the ratio of the rate of parameter changes versus the 
bulk dissipative gap hcq. This separation of time scales, 
naturally provided due to the non-evolving subspace, pre- 
vents the protected decoherence-free subspace from ever 
being left, a phenomenon sometimes referred to as the 
Quantum Zeno effect [33 . 

Equipped with this understanding of the role of the 
dissipative gap, we now demonstrate how to move a Ma- 
jor ana fermion adiabatically along the wire in our dissi- 
pative setup (cf. Fig. |4]a)), which is the key ingredient 
to perform braiding. As an example, we describe the 
move of the right unpaired Major ana fermion from 
site TV {jR = C2n) to site - 1 (7i? = C2N-2)' To 
achieve this purpose, we consider an adiabatic change of 
the last Lindblad operator d^-i of the form: dN-i{0) = 
^[ajy - ttAT + cos6>(ajy_^ + aN-i) - smO{a\^ + ^at)], 
where adiabatically varies with time from ^ = where 
^Ar-i(O) = ttAT-i to = 7r/2 where aAr-i(7r/2) = —aN- 
At the end of this evolution, the site N is empty (vac- 
uum), and the right Major ana fermion moves one site to 
the left. This is the analog of locally tuning the chemical 
potential in the Hamiltonian setting to move Majoranas 
[To]. For the evolution of the Majorana mode popula- 
tion induced by a finite ramping velocity for time T, we 
find {C2N-2C1) = (c2ArCi)exp ^—2 dtO'^/hz^^ describing 
a weak dephasing of the Majorana mode (see appendix). 

Operating this mechanism on a T-junction \T0* in or- 
der to exchange the two modes adiabatically while per- 
manently keeping them sufficiently far apart from each 
other, the unitary braiding matrix describing the pro- 
cess is Bij = exp (f 7i7j) for two Majorana modes z, j, 
demonstrating non-abelian statistics since [Bij^Bjk] 7^ 
for i ^ j. While knowledge of the braiding matrix for 
each pair of Majorana modes is in principle enough to 
demonstrate non-abelian statistics, braiding of a single 
pair of Majorana modes is not physically observable since 
super-selection rules dictate a diagonal density matrix for 
a single complex fermion, which therefore is insensitive 
to the acquired phase. 

We follow ^34* i35j to construct an interferometric ex- 
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FIG. 4. a) Elementary Majorana move: transfer of a Ma- 
jorana edge mode one site to the left by adiabatic change of 
the parameter in the dissipative Liouville operator, b) Illus- 
tration of the initial state p^, braided state p$ (upper row), 
the same states in the changed basis together with the results 
of the fermionic number measurements in the interferomet- 
ric gedankenexperiment. The non-zero correlations between 
Majorana fermions are shown as semitransparent lines. 



periment that explicitly shows the non-abelian nature of 
dissipative braiding of Major anas fermions. The setup, 
cf. Fig. [4] b), is given by two finite wires which host 
four unpaired Majorana modes, 7^,1 and jr^i for the 
left and right unpaired Majorana fermions of the first 
wire, and similarly 71,^2 and jr^2 for the second one. 
These four real fermions correspond to two complex 
fermions aj = {jrj + i7L,j)/2 with j = 1,2. For 



the subspace with an even number of complex fermions, 
we define the following basis |0) = |00) = |vac), and 
|1) = 1 11) = a J a2 1 vac). We now prepare an intial 
state of the system as |^) = (|0) - |1)) /V2 (the cor- 
responding density matrix is = |^)(^|) such that 
(7i?,i7L,2) = (7i?,27L,i) = i- If we now braid the Ma- 
jor anas jR^i and 7l,i, the initial state p^ transforms 
into p^ = Bp^B^, where B = exp (jJl^iJr^i). To dis- 
tinguish the states and p^ in an occupation number 
measurement, we first make a unitary change of basis via 
B12 = exp ( j7i?,i7L,2) such that 



Bi2P^B12 
Bi2P^B12 



|l))((0|-(l|)=p*, 



(9) 



and then measure the number of fermions rii = {a[ai) 
and 712 = (<^2^2) on the wires. In the first case, we get 
with equal probabilities either nj = or = 1, while 
we always get rij = 1 in the second case, cf. Fig. [4]b). 



TOPOLOGICAL ORDER OF THE STEADY 
STATE 



We now show that the robustness of the edge modes 
seen above is indeed related to the existence of topolog- 
ical order in the bulk of a stationary state of Liouvillian 
evolution, and construct a topological invariant, which 
characterizes topologically different states. This classifi- 
cation does not rely on the existence of a Hamiltonian or 
on the purity of a state (in contrast to existing construc- 
tions involving ground states of Hamiltonians) , and can 
be entirely formulated in terms of a density matrix. 

In an infinite system, a stationary state of a Gaus- 
sian translationally invariant Liouvillian is described by 
a density matrix p = n/c>o Pk^ where p/c is a 4 x 4 hermi- 
tian unit trace matrix, which describes the momentum 
mode pair The topologically relevant information 

is encoded in the 2x2 block p2k of the density matrix 
in the subspace with even occupation of the modes 
{a\ak) + {oL^_]^ci-k) = 0, 2 (see appendix). The matrix p2k 
is proportional to ^(1 + n/^a), where a is the vector of 
Pauli matrices and fik is a real three-component vector 
< < 1. The pure states correspond to p\ = pk^ i.e. 
|n/c| = 1 for all /c > 0. We can naturally extend the defi- 
nition of fik to negative k following the change of p2k re- 
sulting from the transformation of the basis vectors in the 



subspace under k 



-k: n^'^ 



x,y 



and 



nt. 



Note that the vector is continuous at A: = 0, ±7r be- 
cause n/c=o and nk=±7r have only z-component. 

Once the vector fik is nonzero for all /c, the normalized 
vector n/e = |n/e|~^n/c defines a mapping ^ 5^ of a 
circle 5*^ (the Brillouin zone — tt < k < tt with identified 
end points k = ±7r due to usual periodicity in the recip- 
rocal lattice) into a unit sphere 5*^ of end points of n/^. 
This mapping, however, is topologically trivial (the cor- 
responding homotopy group 7ri(5'^) = 0), since a circle 
can always be continuously deformed into a point on the 
sphere. We therefore need an additional constraint on fik 
in order to introduce a nontrivial topology. In our set- 
ting, motivated by Kitaev's model Hamiltonian [7], the 
constraint is provided by the chiral symmetry [36l [37] . 
In terms of the density matrix, the chiral symmetry is 
equivalent to the existence of a /c-independent unitary 
matrix S with = 1, which anticommutes with the 
traceless part of the density matrix (n/^a in our case): 
S fik^Tj = —fik^' After representing the matrix S in the 
form S = a<j, where a is a constant unit vector, the chiral 
symmetry condition reads n^a = 0, i.e., the vector fik is 
orthogonal to a for all k. The end point of fik is now 
pinned to a great circle on the sphere such that the 
vector fik defines a mapping from the Brillouin 

zone into a circle. The corresponding homotopy group 
is now nontrivial, 7ri(5'^) = Z, and such mappings are 
divided into different topological classes distinguished by 
an integer topological invariant (winding number). The 
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FIG. 5. Visualization of the topological invariant v for chi- 
rally symmetric mixed states. Left panel: Chiral symmetry 
constrains to a great circle. For pure states, it is further- 
more pinned to unit length (red) . Tuning the Liouville param- 
eters destroys the purity and deforms the circle to an ellipse 
(blue). A phase transition occurs when the ellipse shrinks to 
a line (black). Crossing the transition, the topological invari- 
ant changes sign from +1 to —1. Right panel: Purity \nk{0)\ 
(solid) and average occupation rife (^) = \{l — nz,k{0)) (dashed 
and dotted) for various fixed transition parameter for both 
sides of the transition (see text). Throughout the transition, 
the system is half filled. 



explicit form of this invariant reads (see appendix) 



^ = — / dka- {fikX dkfik) e Z. (10) 

Geometrically, u counts the number of times the unit vec- 
tor n/e winds around the origin when k goes across the 
Brillouin zone, cf. Fig. [5] Importantly, the winding num- 
ber u distinguishes topologically different density matri- 
ces for translationally invariant Gaussian systems with 
chiral symmetry without restriction on the purity of the 
state. 

Let us now discuss specific examples of the steady 
states resulting from quasi-local dissipative dynamics. As 
shown in the appendix, for momentum space Lindblad 
operators jk = C^^/c, = (^fc, v/c), = (^/c, with 
(unnormalized) Bogoliubov functions Uk^ Vk which do not 
induce a current in the steady state ((a^a/c) = (a^^a_/e)), 
the solution for the vector fik reads 



1 



i^k^^l^ik^ (11) 



2 1 ■■-y.'^ ■■''y,-K 
rriz^k ^rriz-k 



with hZk = {(.k^k + ^!_ze^^-/e)/2. For quasi-canonical de- 
formations (cf. Eq. ([6|), the additional symmetry prop- 
erty U-k = and V-k = —"^k simplifies the solution 
to fik = m/c, implying the purity of the steady state, 
\nk\ = 1 for all k. The solution corresponds to the ground 



state of some Hamiltonian. For = Og = 7ts/2 with an 
integer s, the system has a damping gap closing point, 
and the vector fik has only z-component (—1 for even s 
and +1 for odd s) for all /c, leading to u = 0. For other 
values of ^, one has u = ^1. However, as one cannot de- 
fine the direction of a for ^5 , a global definition of u for all 
6 is not possible, with the consequence that the potential 
change of the sign of u when passes Os is meaningless. 
Note that the gap closing points do not correspond to 
a phase transition here, because one has either a com- 
pletely empty or completely filled lattice, such that no 
thermodynamic observables of a phase transition could 
sensibly be defined. 

For non-canonical deformations, Eq. ([7|), visualized in 
Fig. [sj the steady state density matrix is mixed, \fik\ < 1 
(cf. also right panel in Fig. [5|. Importantly, the topo- 
logical order persists for quite strongly mixed states, and 
we find again u = ±1 for ^ Og- The difference to the 
previous example is that at 6> = 6>s , not only the direction 
of a but also the topological invariant is not defined: fik, 
alined in the ^/-direction for all /c, has zeroes: nfe=o,7r = 0, 
meaning physically that these modes are in a completely 
mixed state. The "loss" of topology at ^ = 6^5 can be 
viewed as a non-equilibrium topological phase transition 
[131 El EH] as a result of changing the Liouville parame- 
ters: The system has well defined thermodynamic prop- 
erties, since it is half filled, fi{0) = J ^fikiO) = 1/2 for all 
0. The closing of the dissipative gap at Og leads to critical 
behavior, which manifests itself via diverging time scales, 
resulting e.g in an algebraic approach to steady state (as 
opposed to exponential behavior away from criticality) 
[2Q[ [2T| [23l [39] . Importantly, the vector fik in the steady 
state has the reflection property fik{—S6) = Pfik{^SO) 
for ah k, where P = diag(l, 1, -1) and 60 = 0-Os. There- 
fore, the symmetry pattern of the steady state is identical 
on both sides of the transition, ruling out a conventional 
Landau-Ginzburg type transition and underpinning the 
topological nature of the transition. 

We see that in our system the stationary state of Liou- 
villian evolution in an infinite system is indeed character- 
ized (for ^ Os) by a nontrivial topological order. If the 
system is finite, the edges separate the topologically non- 
trivial bulk from a vacuum, which is topologically trivial. 
Similar to the Hamiltonian case of topological insulators 
or superconductors [6 , the "jump" in the topological or- 
der guarantees the existence of edge states. For a pure 
stationary state (quasi-canonical case) this follows from 
the formal analogy with the Hamiltonian case in which it 
is well-established (see e.g. Ref. [40 ). More generally, for 
a mixed state (non-canonical case), the arguments follow 
the line of Ref. [7 using an alternative equivalent form 
of Eq. (10) for the topological invariant discussed in the 



appendix, cf. Eq. (41 ). The topological origin of the edge 



modes explains their stability demonstrated above. 
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IMPLEMENTATION IN COLD ATOMIC GASES 

Here we devise an interacting (quartic in the fermion 
operators) Liouville operator which at late times reduces 
to the quadratic Major ana Liouvihe operator Eq. ([2]). 
The late time limit and the steady state here play the 
role of a low energy limit and the ground state in equi- 
librium, respectively, where the physics of weakly corre- 
lated superconductors is universally described in terms of 
Bogoliubov quasiparticle excitations for a variety of mi- 
croscopic models. Our prescription may thus be seen as 
a microscopic "parent Liouvillian" , providing one possi- 
ble microscopic realization of a Liouville operator which 
generates the desired properties at and close to steady 
state. 

The implementation idea follows closely an earlier pro- 
posal for bosons [20 and is based on a conspiracy of laser 
driving and engineered dissipation made possible by im- 
mersion of the fermionic target system into a superfluid 
BEG reservoir. The interaction with the bosonic dis- 
sipative reservoir is microscopically based on a conven- 
tional s-wave fermion-boson density-density interaction 
and stands in marked contrast to the proximity effect to a 
BCS superconductor exploited in solid state implementa- 
tion proposals of the (Hamiltonian) Majorana wire. Our 
scheme is illustrated in Fig. |2] and explained in more 
detail in the appendix. It yields the following number 
conserving Liouville dynamics. 



(12) 



These Lindblad operators give rise to dissipative pairing 
in the absence of any conservative forces, a mechanism 
based on an interplay of phase locking and Pauli blocking 
established recently |22J. The relation to the Majorana 
operators is apparent in the thermodynamic limit, where 
it can be shown that the following general relation be- 
tween fixed number (J^) and fixed phase {ji) Lindblad 
operators holds (see appendix). 



(13) 



Here Cj{Ai) are creation (annihilation) parts, respec- 
tively, which for Eq. (12) read Cj = ^(^J + aj^-^), = 
|(ai — a^+i). In consequence, we see that indeed pre- 
cisely Kitaev's quasiparticle operators are obtained as 
Lindblad operators, ji = ai. The description in terms of 
fixed phase operators becomes appropriate at late times, 
i.e., close to the steady state, where an ordering princi- 
ple is provided by the macroscopic occupation of only a 
few correlation functions. The explicit calculation shows 
that Eq. ([2| is produced with effective dissipative rate 
Hi = k/S (see appendix). There we also discuss the lead- 
ing imperfections, showing that they preserve the chiral 



symmetry necessary to remain in the above described 
topological class. 



CONCLUSIONS 

In this work, we established a complete list of topolog- 
ical features familiar from ground state physics of cer- 
tain Hamiltonians in engineered dissipative dynamics, 
highlighting the universality of the concept of topolog- 
ical order in quantum mechanical many-body systems. 
While the present work has focused on the conceptually 
simplest system of a quantum wire, the idea of dissipa- 
tively induced topological order is more general and will 
be present in other physical systems and in dimensions 
higher than one. 
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Appendix 
Implementation in cold atomic gases 

Microscopic Model - Here we specify a physical setup 



leading to Eq. (12), following an earlier proposal for 
bosons 23 illustrated in Fig. [5] We start from an 
optical superlattice setting, with lower sites which make 
up the fermion wire and correspond to the lowest Bloch 
band of the lattice, and auxiliary sites associated to the 
second Bloch band located on the links. These two Bloch 
bands are coupled via driving lasers, whose Rabi frequen- 
cies ±Q are chosen with opposite sign for each pair of 
lower sites. This amounts to a commensur ability condi- 
tion of driving and lattice laser, and ensures the relative 
minus sign in the annihilation part of the Lindblad op- 
erators of Eq. (12), Ai oc ai — a^+i. By immersing the 



whole setting into a BEG reservoir, particle superposi- 
tions in the upper band can spontaneously decay back 
to the lower one by emission of a Bogoliubov phonon 
into the BEG bath. This process is isotropic and short- 
ranged for suitable bath parameters, giving rise to the 



two-site creation part in Eq. ( 12 ) with relative plus sign, 
Ci (X ai -\- a^+i. Microscopically, this dissipative mecha- 
nism requires a standard s-wave fermion-boson interac- 
tion between system (fermions in the optical superlattice) 
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and bath (BEC) particles [20]. An effective single par- 
ticle microscopic model such as Eq. (12) is obtained by 
integrating out the upper band under suitable detuning 
conditions for the driving laser. Using the recently devel- 
oped tools of single site addressability for optical lattices 
[26| l27] , an edge can be constructed by cutting the su- 
perlattice in the right place cf. Fig. |2] 



Eq. (12) describes an interacting, number conserv- 



steady state properties in the late time dynamics. The 
resulting mean field dynamics close to the steady state is 
given by 

dtp = ^i^k {jkpjl - HjkJk^ p}) , (14) 

k 



ing Liouville dynamics, which leads to dissipative pair- 
ing [22], described by a pure BCS-type wavefunction for 
paired fermions. There are two ways of seeing this, either 
by explicit construction of the dark state wavefunction 
(fixed particle number) or by deriving a suitable mean 
field theory (fixed phase). In the thermodynamic limit 
both procedures are equivalent. Here we follow the sec- 
ond option, and highlight the more precise connection 
between the two approaches below. 

Mean field theory at late times - At late times, follow- 
ing Ref. [22 we can make use of the proximity to the 
steady state to derive a quadratic mean field theory for 
this dissipative dynamics. In the spirit of BCS theory, 
the ordering principle behind this approximation is the 
macroscopic occupation of only a few correlation func- 
tions in steady state, which can be evaluated explicitly 
on the exactly known steady state. To implement the ap- 
proximation, as usual in BCS type mean field theories, 
we give up exact particle number conservation and work 
with a fixed phase. 

We thus start from number conserving Lindblad op- 
erators Ji = CjAi, where the creation part Cj = 



where jk are fermionic quasiparticle operators up to a 
normalization, obeying the anti-commutation relations 

{jkJl} = + \v,\^)Sk„{jk,j,} = {jijl} = 0. In 

our example, \uq\'^ + \vq\'^ = 1 and = i.e. the dy- 
namics has constant damping rate for all modes. In fact, 
transforming back to position space produces precisely 



= 0, we see that the 



Eq. ^J^ = Ek^~''''jk=U^^ 

Since jk{uk + e'^Vka'^_^al)\Ydic) 
above ansatz indeed provides the correct steady state 
solution. In addition, the equivalence of fixed number 
and fixed phase wavefunctions can now be justified in 
our nonequilibrium context a posteriori in the thermo- 
dynamic limit: the fixed phase BCS state has relative 
number fiuctuations AA^^ — )-(^) 



where TV 



is the particle number operator and N the number of 
degrees of freedom. 

We emphasize that, remarkably, the late time evolu- 



tion of the master equation (12) naturally gives rise to 



and the annihilation part Ai 



with translation invariant complex position space func- h.c. in Eq. ^\A). 
tions Vi-j^Ui-j] in the example Eq. (|3| discussed in 
the text, Vi-j = \{5ij + 5ij^ij),Ui-j = ^{Sij - 
Si-^ij). For the practical calculation, we switch to 
momentum space, where we have Jk — ^^e^^^^Ji = 
^qC]^_kAq. The Fourier transforms for the creation 
and annihilation part are local in momentum space, 
Vka\,Ak = Y.^e^^'^'Ai = Ukau, 



quasilocal squeezing of fermions. The mechanism behind 
this effect is in close analogy to a superconductor: the 
system here acts as its own reservoir by providing an or- 
der parameter, which allows for the appearance of the 
off-diagonal pair annihilation / creation terms ~ aka-k, 



General relation of fixed number and fixed phase Lind- 
blad operators - Here we show Eq. (13). Suppose we 



are given a fixed phase Liouvillian defined by a set of 
Lindblad operators of a form 



ji = cj+A 



(15) 



' COS ■ 



ie ^^/^ sin | in the example 



with Vk 

(due to the structure of the Liouville operator, the expo- 
nential prefactors are irrelevant and can be omitted). In 
the following we work with functions with the property 
U-k = ±.Uk^V-k = T'^/c, which however do not necessar- 
ily obey a normalization constraint. 

For a setting with fixed phase, we now make an ansatz 
for the steady state density matrix of the product form 
p = Y[k>oPk^ with density matrix for each mode pair 

Pk = {uk^e'^Vka'if^al)\Ydic){vdic\{ul^e~'^vlaka-k)' In- 
serting this ansatz into Eq. ([l2| in momentum space, 
and using the projection prescription pk = tT^±kP on 
the mode pair we obtain the equations of motion 
for the single pair density matrices in the presence of 
nonzero mean fields. These result from the coupling to 
other momentum modes, with values determined by the 



using the above conventions, such that the momentum 
space Lindblad operators read jk = u^ak + Vko}_k' The 
dark state (with property ji\BCS,0) = jk\BCS,0) = 
for all i or /c) is conveniently formulated in momentum 
space and reads 

\BCS,0) (xexp(e^^G^)lvac) = ]J(1 + e^V/c^i-fe4)|vac), 
= ^ ^ (fikalual, (fik = — = -^-k- (16) 



k 



Uk 



Now we want to show that the number conserving version 
is given by the following expressions for the Lindblad 
operators and related dark state wavefunction. 



Ji = CjAi, \BCS,N) (X G^^lvac). 



(17) 



To prove it, we proceed in momentum space. For the 
normal ordered Lindblad operators J^ = ^qC^-k^gi 



9 



the dark state property Ji\BCS,N) = Jk\BCS,N) = 
for all i or is equivalent to the vanishing of the following 
commutator for all /c, 

[J/e, G"^] = Vq-kUq(Pqal_^alq (18) 

q 

This is true if and only if ^^^^^ = for ah k, i.e. for 

^ UqVq-k ^q-k ' 

the wavefunction (pq = Vq/uq as claimed above. 

Note that working with the functions \BCS,N) for 
finite N requires an infrared momentum cutoff ql ~ l/L 
such that N = nL {n the density of particles), that is 
consistent with the thermodynamic limit A/" — ^ oo,L ^ 
oo, n = N/L const. 

Thus, for a given real space quadratic master equa- 
tion, we can immediately construct a number conserving 
version and vice versa by the above arguments, and in- 
dicate the respective fixed number of fixed phase exact 
dark state wavefunctions. 

Imperfections - In the above implementation scheme, 
the annihilation part Ai = ai — a^+i is well under con- 
trol since it relies on a locking of lattice and driving 
laser. The dominant imperfection appears in the cre- 
ation part: In the case that the phonon wavelength in 
the BEC bath is not smaller than the lattice spacing 
(subradiant case), the decay may take place over sev- 
eral lattice sites, such that e.g. Cj = ^(aj + al_^-^ + 
e(a|_-L + <^i+2)' giving rise to fixed phase Lindblad oper- 
ators jk = ""^^^[(cos I + ecos ^)a^^ + i sin |a/e] with 
Nk = 1 + 2e cos I cos ^ + cos^ ^ in momentum space, 
which exhibit fermionic anticommutation relations guar- 
anteeing a pure steady state. The resulting vector fik 
still lies in a plane for arbitrary e, i.e. the imperfection 
preserves the chiral symmetry. 

Properties of the finite size system 

Equations of motion in the Majorana basis - For the 
practical treatment of the finite size system quadratic in 
the fermion operators and with Gaussian initial condi- 
tions, it is convenient to encode the information in the 
covariance matrix of second moments in the real Majo- 
rana basis. For convenience, we repeat and extend here 
the definitions of the main text. The Majorana opera- 
tors obey C2j = (a] + a^), C2j-i = i(a] - a^), c] = cj, 
{ci^Cj} = 2Sij^ where j = 1, ...A/" labeling the physical 
sites in a one dimensional lattice. Following [23], we start 
with a quadratic master equation written in the Majo- 
rana basis, i.e. dtp = -i[H,p] ^ i^Y^iUpjl - \{j\ji^p} 
with 1-L = He for 2N x 27V hermitian matrix H and 
ji = ljc^j\ = c^l*^ for 2N component column vectors 
li^c. The Liouvillian parameters are then encoded in a 



hermitian 2N x 2N matrix M = J2ih ^ ^J- The equa- 
tion of motion for the covariance matrix with entries 
Fab = i([ca,C5]) rcads {hi = 1) 

S,F = = -i[i7,F]-{X,F}-r, (19) 

with real matrices X = 2ReM = , Y = 4ImM = 
—Y^. The spectrum of X is positive semidefinite [24] 
and a pure state is signaled by the eigenvalues of F^ being 
all equal to —1. 

Zero modes - The real symmetric matrix X is diago- 
nalizable and has real eigenvalues and -vectors. First we 
show the existence of zero eigenvalues for this matrix for 
a wide class of finite systems via explicit construction of 
the corresponding eigenvectors. The ideal Lindblad op- 
erators for the dissipative topological ID quantum wire 
are characterized by the following vector in the Majorana 
basis: 

^J=(0, 0, (/,)2„ (/,)2,+i, 0, 0), (20) 

with {lj)2j = i and (/j)2j+i = 1- Let us now consider the 
more general vector Ij specified with 

i}j)2j-i = e2j-i, {lj)2j = i(l + e2j), ^^^^ 
(^j)2j+l = (1 + e2j + l), {lj)2j+2 = ie2j+2. 

For each /j, there are two orthogonal vectors given by 

; X (22) 

vIj={0, 0, 0, 1, O), 

and by construction, they fulfill ij • vlj = ■ VL,j = 
ij ' ^Rj = • vrj = 0. Due to the structure of the 
matrix M as a sum of outer products of the Ij vectors, 
we find precisely two zero modes for any finite system 
size N given by the 2N component vectors 

N 

7L = X]^L,jC2^-l, 

J = l 

N 

IR = X]^i?,jC2j, 
J = l 

yRJy^2N = Yi \~T~rT^^ ) ' ^RJ=2N = 1, 

(23) 

where Jl^r represent the left / right Majorana modes. 
The reason behind these two vectors being zero vectors 
of the matrices M and M* is that they are orthogonal 
to the {2N — 2) linearly independent vectors Ij and /* in 
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a space of dimension 2N . Hence, vl^ vr^ {Ij} and {/*} 
form a complete basis in this space. 

Let us now focus on the deformations studied in the 
main text. Here, the Major ana modes are characterized 
by vectors 

(i) Canonical deformation: 



7 — 1 N—j 



(24) 



(25) 



(ii) Non-canonical deformation: 
with e = sm ^-cos^ ^ rpj^ associated localization length 

cos y+sm ^ 

is then given by l\oc = —ci{\og\e\)~^ ^ with a the lattice 
constant. 

In the case of a deformation being non-homogenous 
but with value changing randomly from site to site, the 
zero modes are still present, but the expression for 
and ^R is given by 

3 3 
'^L,2j-1 = fjci, VR^2(N-j+l) = (26) 
i=l i=l 

Dissipative isolation of the subspace - As argued in the 
main text, in addition to the existence of a zero mode 
subspace, the dissipative evolution must ensure its iso- 
lation from the bulk. The conditions for such a situa- 
tion are readily formulated generally, without reference 
to the quadratic setting. For a set of Lindblad operators 
Ji making up the total Liouvillian C[p] = Jipjj — 
we may introduce projectors on the edge 
(zero mode) and bulk subspaces, p and q = 1—p^ respec- 
tively. A decoupled edge subspace appears if the Lind- 
blad operators Ji are block diagonal, Ji^pq = Ji^qp = 0, 
with the edge block identical to zero, Ji^pp = 0. We then 
obtain a dissipative evolution for the density matrix in 
this projection. 



dt 



ppp ppq ^ _ 
Pqp Pqq 



1 t 

~2Pp(l'^3,qq'^^: 



(27) 

The bulk dissipative evolution jCj^qq[pqq] = Jj,qqPqqj]^qq — 

^{'^j qq'^jyqq^ pqq} Lindblad form. The density matrix 
in the p subspace ppp is a constant of motion. The cou- 
pling density matrix elements pgp, h.c. damp out accord- 
ing to Pqp = e~^o '^Iqq^j'«'?^p^p(t = 0), i.e. exponentially 
fast in the presence of a dissipative gap. 

This situation is indeed present in our quadratic set- 
ting. We switch to the spectral representation X = 
^r\r){r\^ where ^ are the nonvanishing eigenval- 
ues and |r) the corresponding eigenvectors. In addition, 
we have zero modes Aq, = which do not contribute to 
the spectral decomposition, with eigenvectors \a). With- 
out loss of generality we choose orthonormal eigenvectors 



{a\b) = 5ab' In this basis, the equations of motion read 



dt 



^a(3 Lq/c 







-(Xr 



as 



Y 



(28) 



where Yab = {a\Y\b). Following the above discussion, the 
zero modes of the matrix X are also zero modes of the 
matrix Y and vice versa. This implies Y^s = = and 
shows the decoupling of the edge and bulk subspaces, as 
well as Yai3 = 0, which defines the zero modes subspace. 
This reflects the structure of Eq. (27) for a quadratic 
theory in the Majorana basis. 

Adiabatic parameter changes - We study the evolution 
for a time dependent Liouville (or Hamilton) operator in 
the instantaneous basis in the quadratic setting, where 



the form of ( 28 ) is kept but the eigenvalues and -vectors 



of X now depend on time. The transformation into the 
instantaneous basis \a{t)) = U{t)\ao), where |ao) is the 
initial reference basis, is now orthogonal. The connection 
is given by the matrix A = {a{t)\b{t)) = {ao\U^U\bo)^ 
which is antisymmetric due to normalization, dt{a\b) = 0. 
With Hamiltonian in the instantaneous basis given by 
hab = cind hermitian matrix A = iA, the equation 

of motion in this basis reads 



dtT = -i[h^A,T]-{X,r}-Y. 



(29) 



As already mentioned in the main text, the shift in the 
Hamiltonian h ^ h -\- A shows the emergence of a gauge 
structure due to the explicit time dependence of the 
eigenbasis, irrespective to whether the time evolution of 
r is generated by a Hamiltonian or Liouvillian. 

We study the simple example of local adiabatic pa- 
rameter changes from the main text: First, we consider 
a system of two physical sites, specified by the vector 
k = ^(0, cos^, i, — sin^). The spectrum of X is then 
doubly degenerate with eigenvalues 0,1/2. The connec- 
tion matrix is ^4 = 0{l — az)^o-y. The resulting equations 
of motion for the elements with nonzero rhs are, 



(30) 



dt^2?> — — ^ + ^r34, 9tr34 = — — ^r23 — 1. 

The system is easily solved using adiabatic elimination 
in the limit where the parameter changes are adiabatic, 
i.e. d(t) <C 1 for all times. In particular, we then obtain 
for the slowly evolving variable in the zero eigenvalue 
subspace 
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-2rr 



12 



(31) 



with solution (written with dimensions restored) 

ri2(t) = ri2(0)exp I -Ik-' \ dt'O' ) , (32) 
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describing a dephasing of the Major ana mode. Actually, 
for linear time changes 6{t) = at the integral approaches 
o^T/k <C 1, i.e. the Majorana mode population is only 
weakly affected. 

Now, we observe that this result also applies to the 
general case of N sites. This is due to the structure of 
the steady state resulting from the ideal case described 
by Eq. (|3| with Majorana modes paired between different 
physical sites. The only time dependent changes occur on 
the Majorana sites 1, 2N -2, 2N -I, 2N , while the rest of 
the system remains decoupled (the situation is visualized 
e.g. in Fig. [T]). Therefore, the result (32) is general. 

We thus show that the Majorana subspace remains 
protected for adiabatic parameter changes also in the dis- 
sipative case. The argumentation will be valid and useful 
in any situation where there is a gap in the damping spec- 
trum. 



Topological Invariant and Master Equation 

Here we present some details on the master equa- 
tion and topological invariant for a translation invari- 
ant quadratic Liouville operator with chiral symmetry in 
terms of the density matrix. 

Density Matrix - We consider a translation invariant 
setting, where this property applies to both the initial 
state and the Lindblad operators defining the Liouvillian. 
For a quadratic Liouvillian, the density matrix then takes 
a product form in momentum space p = n/c>o P^^ where 
pk is the density matrix of the momentum mode pair 
Using fermion superselection rules, the 4x4 matrix pk 
can be written in a block-diagonal form 



Pk 



Pik 

P2kJ' 



(33) 



where pik = diag(p^^, p^].) is a diagonal 2x2 matrix 

with non-negative elements in the subspace with 
odd total occupation of the modes -\-k and —k equal 1, 
(a|.a/c) + (a^^a_/c) = 1, and p2k is a 2x2 hermitian matrix 
in the even occupation subspace (a^a/e) + (a^^a_/e) = 0, 2. 

Taking normalization (p^^ + + trp2fe = 1) and her- 
miticity into account, the matrix p2k can be written as 



P2k = trp2/c^(l + Q/c), 



(34) 



where Q/c is a traceless hermitian matrix, = fik^ with 
(J being the vector of Pauli matrices and fik a real 3- 
component vector < |n/c| < 1. For a pure state one has 
\nk\ = I for all k. In this case p±k = 0, tTp2k = 1, and p2k 
takes the form of a projector, p2/c = P^k- The situation 
n/e = for all k signals a completely mixed state. 

It is convenient to introduce the 2x2 matrix p2k = 
{tTP2k)~^P2k, such that 



which carries all topological information encoded in the 
density matrix pj.. This is because the other quantities 
entering^,, n.mely p% p% and tvp,„ are just non- 
negative numbers subject to the constraint p^^l + p^_}l + 
trp2/c = 1, and the space defined by these conditions is 
topologically trivial. 

The matrix pj. and, therefore, the vector fik are defined 
for /c > 0. It is convenient, however, to extend them 
also to negative values of k. This can easily be done 
by noticing that under the change k ^ —k the elements 
|0/c, 0_fe) and \ lk^l-k) of the basis in the even occupation 
subspace change as |Ofe,0_/c) |0_/c,0/c) = |0/c,0_fc) and 
|l/c,l-fe) |l-fe,l/c) |lfe,l-/c). As a result, the 

off-diagonal elements of change their sign, while the 
diagonal ones remain unchanged. That is, we have the 
transformations 



O'zQkO'z, 



n-k 



^z'^k^ 



(36) 



with Sz = diag(— 1, — 1, 1). Note that the thus defined 
vector n/e is^ continuous at /c = and k = ±7r because 
Qk=o and Qk=±T^ have only diagonal elements {fik has 
only z-component) due to the fermionic nature of the 
particles. If in addition |n/e| ^ for all /c, the vector 
fik defines a mapping of the Brillouin zone (topologically 
equivalent to S^) into a sphere 5*^, k ^ fik = Uk/ \^k\' 

Topological Invariant - The mappings constructed 
above are topologically trivial because the homotopy 
group TTi of a sphere S'^ is trivial, 7ri(5'^) = 0. In other 
words, every mapping ^ S'^ can be continuously de- 
formed into a trivial one that maps the entire Brillouin 
zone into a point on a sphere {fik = n is /c-independent). 
As a result, a general density matrix of the form (33) 



does not have topological order. The situation is differ- 
ent for the density matrices obeying the chiral symmetry 
[36l|37]. In this case, there exists a unitary matrix H that 
anticommutes with Qk for all /c. 



(37) 



The condition of the chiral symmetry in our case has a 
simple geometrical interpretation. Namely, if we write 
the matrix S in the form S = aa, where a is a real unit 
vector, \d\ = 1, then Eq. (37) is equivalent to 



afik = 0, 



(38) 



i.e. the vector fik for all k belongs to the plane that 
is orthogonal to the vector a. As a result, the vector 
fik determines a mapping of the Brillouin zone {S^) into 
a circle S^^ which is an interception of the unit sphere 
S'^ with the plane passing through its origin (one of the 
great circles). The mappings of this type can be divided 
into topologically different classes that are characterized 



P2k 



^i + Qfc), 



by an integer winding number (10) or (42) because the 



(35) homotopy group in this case is nontrivial, 7ri(5'^) = Z. 
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To construct the topological invariant that distin- 
guishes topologically different density matrices we use 
the fact that Ql = n|l and introduce the normalized 

matrix Qk = \ fik\~ ' 
matrix can always be written as 



Qk with a unit square, Qi = 1. This 



Qk = U^ 







Qk 



U 



-mz 



(39) 



(40) 



1, with 



where 0/^, mx,k^ and ruz^k are real, 
some unitary matrices U and V (in the first case S = 
±U^azU, in the second case S = iVV^V). Therefore, 
all information about the stationary state is encoded in 
the phases (j)k or in the unit vectors rhk = {jrix^k-, 0^ 
in the {x — z)-plane. As a result, the topologically differ- 
ent classes of stationary states with the chiral symmetry 
can be classified by an integer winding number 



W[p\ 



d<pk 



(41) 



dk{m 



-I 

Jbz 



dk 
2^ 



drrix 
~dk 
drrix 
dk 



m. 



druz 
'~dk 
dm 



dk 



where in the last line we extend the integration over the 
entire Brillouin zone k G [— tt, tt] using the continuation 
from positive to negative k as discussed above, cf. Eq. 



(36). We see that W[p] indeed indicates the number of 
times the end of the vector rfik winds around the origin 
when k goes over the Brillouin zone k G [— 7r,7r]. 

In terms of the matrix Q/c, the invariant can be written 

as 



W[p] 



/ dktii^QkdkQk)^ (42) 



Inserting S = aa^Qk = nk^ into Eq. (42) yields Eq. 



(10) in the main text. 



Note that the definition of the topological invariant de- 
pends on the vector a, which can only be defined up a sign 



(S and — S are equivalent). While k independent, this 
vector can in general depend on the parameters of the Li- 
ouvillian. As a result, the global definition of v requires 
the existence of a continuous choice of one of the two 
branches of a in the entire parameter space. Whether it 
is possible or not, depends on the parameter space itself: 
One has to be able to connect any two points by some 
deformation path, along which one can define a continu- 
ously. If this is not the case, the global choice of a and, 
hence, the global definition of v is not possible, and one 
has arbitrarily chosen a in different " disconnected" parts 
of the parameter space. The relative sign of v in different 
parts has then no significance. In any case, however, the 
nonzero value of v indicates nontrivial topological order, 
irrespective to its sign. 

Equation of Motion and Solution - Using the product 
form of the density matrix, applying the projection pk = 
tr^±/eP (all mode pairs but ±/c are traced over) we obtain 
the equation of motion for each pair, which reads {n = 1) 



dtpk = ^^kpkjik - ibikjcTk^pk}^ 



(43) 



cr = ± 



jk UkCk + Vkci^ ik^k, & = 



Uk 
Vk 



Ck 

J 



for k ^ 0, ±7r. Here Uk^Vk are arbitrary complex func- 
tions of momentum; in particular, they need not to fulfill 
a normalization condition. For the solution, we study the 
set of single particle correlation functions (covariance ma- 



trix) which is closed due to the quadratic nature of (43). 
Using suitable symmetrizations, the equation of motion 
for the correlation functions can be written in terms of 
four component vectors 



dtNk = -Kk{{l + Ak)Nk - M|), Ak 







(44) 



where the vector Nk is defined in terms of the single par- 
ticle correlation functions, and we use the further defini- 
tions 
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Mt 



l^k 



^0,k 



ml 



^li,k 



^{Ck — C-k) \ 

2 (Cfe + C-k) 



rrix^k 

^y,k 



J 



-k 



Ck = {ckc\) - {c\ck), 



/i = 0, ...,3, Kk 



2 



y,-k 

\ rriz^k ^rriz-k ) 



m 



'0,k 



ml 







- rho-k), 



rrix^k 



-m^ 
- m,, 



y,-k 

\ mz^k - rriz-k J 



(45) 



and we used (1) = 1. We summarize some properties of 
these equations: 

(i) By construction, all components of iV^, M^^^ are even 
or odd eigenfunctions under momentum reflection. In 
particular, 

n_/e Szfik, m^-k Szml, rh% = -Szml. (46) 

(ii) Despite the appearance of four components in Eq. 
( [44| ), all physical information on the steady state can be 
stored in the vector n^. This is because given this knowl- 
edge, the first line in this equation only carries redundant 
information: no,/e = /c ~ ^k^k- "^^^ general time de- 
pendent solution may need four real numbers for each 
mode pair ±/c for a complete description. 

(iii) Fermionic operators - Important differences between 
the steady states occur depending on whether the Lind- 
blad operators describe fermionic quasiparticle operators 
(up to a normalization). More precisely, the anticommu- 
tation relations are 

{jk.jl}=iUkSk,,, (47) 

{jkjq} = {UkV-k^U-kVk)h-q, h. C. 

Introducing jk = jk/ {^l^kY^'^ ^ the first anticommutator 
can be normalized if ^^^^ is strictly positive. The Dirac 
algebra is then fulfilled if the equal charge operators an- 
ticommute, i.e. for 

UkV-k = -U-kVk' (48) 

If these conditions are fulfilled, we call the deforma- 
tion from the ideal Lindblad operators Eq. ([3| quasi- 
canonical^ as it is related to a true canonical transforma- 
tion by a simple rescaling. (In the specific quasi-canonical 
deformation considered in the main text, Eq. (|6]), even 
the stronger conditions U-k = —Uk^^-k = ^k apply.) In 
this case, the solution simplifies since the components 
m^^/c themselves are odd or even eigenfunctions under 
momentum reflection, notably m_/e = Szfrik. We then 
have Mfe^ = {0,ml),M^^ (0,0^). In this case, the 
steady state solution simplifies to Hk = rhk = ^^^^/c/^o,/c 
fulfilling the purity condition |n/c| = 1 for all /c, such that 
the density matrix takes the form of a projector, p\ = p^. 



This density matrix then describes the ground state of a 
Hamiltonian H = ^j^H^^Hk = ^^^mfca^fc. In con- 
trast, when the vector rhk encoding the structure of the 
Liouvillian does not have the above transformation be- 
havior, m-k 7^ SzTUk for some modes, the steady state is 
mixed. 

(iv) Current-free situation - For no,/c = the sys- 
tem carries not current, {Jk=o) = 0, where the 
current operator in momentum space reads Jk = 
e~^^/^ sing'a^_^^2^g+fc/2- Solutions of this type are 

obtained for the cases Mlj. = or mo,-fe = ^o,/c5 

and M-k ■ Mj^ = 2mz^-k'^z,k- The first condition has 
a simple interpretation in terms of symmetric damping 
of ±/c modes, similar to the symmetry in the excitation 
spectrum of a Hamiltonian in the absence of a magnetic 
field. In this case, the steady state solution simplifies 
to n/e = m|,mQ = 0, and in general describes a mixed 
state |n/e| < 1, which cannot be written as a ground state 
of a Hamiltonian. We note that for fermionic jk (up to 
normalization), this condition is always fulfilled and the 
system is current-free. 

(v) The solution of the linear equation is given by 

Nk{t) = e-^^^7Vfc(0) + (1 - e-^''^^')L-^Mk (49) 

with Lk = /^/c(l -\- Ak). Clearly, symmetries of the time 
evolving density matrix are determined both by the ini- 
tial state A^/c(0) and the Liouville operator parameters, 
i.e. Lk^Mk- For a unique steady state (all eigenvalues 
of L/C strictly positive), all memory of the initial state is 
lost and thus the symmetry of the Liouvillian is inher- 
ited by the steady state density matrix. If a symmetry is 
shared by both initial state and Liouville operator, that 
property is conserved during the evolution. In particular, 
we note that generic initial states like Gaussian thermal 
states are chirally symmetric. 

(vi) The damping spectrum for the correlation functions 
is determined by the eigenvalues of />:/c(l + A^). In the 
case of the Lindblad operators being fermionic quasipar- 
ticle operators, it is fourfold degenerate and coincides 
with tik- More generally, the spectrum is given by 

\]:^ = Kk, = Kk{l±\ml\), (50) 
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with the first eigenvalue doubly degenerate. Positive 
semidefiniteness is seen from |m^| < 1. 

Examples - In the main text we discuss two examples 
with quasilocal (nearest neighbour) Lindblad operators, 
the quasi-canonical and the non-canonical deformation. 
In momentum space, their extended versions are given 



by jj^^ = ^fc^^'^/c with £ = c^n for the canonical and non- 
canonical deformation. 



V2 



1 

V2 



-ie^^ siiiO sin | 
cos 9 cos I 



(51) 



e^^(cos l9e-^^/2 _ sin(9e^^/2) \ 
cos 6>e^^/2 + sin l9e-^^/2 J ' 



steady state solutions 



-(c) 



An) _ 



1 



1 + cos(2^) cos k 

— sin (j) sin k 
cos (j) sin k 
- sin (2(9) cos (k)j 



sin sin (26) sin/c 
cos ^ sin {20) sin k 
— (cos {29) + cos k) y 



(52) 



and the damping spectra for the correlation functions are 
specified with (/^ = 1; cf. Eq. (ISOl)) 



(c) 



1 + COS {29) cos /c. 



.o(c) I 



0, 



(53) 



I -^Ci{n) I 



= I cos(2^) cos 



Both cases reduce to the ideal case Eq. (|3| for (9 = 
7r/4 + STT, (/) = (5 integer), and feature damping gap 
closing points at 9 = si:/ 2. The steady state in both 
cases is current free. The nature of these gap closing 
points is conveniently discussed by studying the rela- 
tions of the steady vector fik in their vicinity. In par- 
ticular, we find the transformation rules (with parame- 
terization 9 = 7t/2 -\- S9^ and restricting to the interval 
7r/4 < 9 < 37r/4 without loss of generality) 

for arbitrary (j) with Sz = diag(— 1, — 1, 1). Thus, the 
steady state vectors on both sides of the transition point 
relate by reflections of two of their components with the 
third one fixed for all modes /c, leading to the conclusions 
drawn in the text. In addition, from the transformation 
properties alone we can deduce basic statements on ther- 
modynamic properties of the system at the transition 
point. In the first case, the situation is thermodynam- 
ically trivial because the system is either fully filled or 
empty: at the transition point 69 = 0^ hy Eq. (54) the 
vector points in ±z-direction for all k. The additional 
constraint of purity only allows for fik = (0,0, ±1) for 
each k. A steplike change of the sign for some values of k 
(Fermi surface) is not possible in our context, since such 
step function in momentum space could only be synthe- 
sized from Uk^Vk which are highly nonlocal in position 



space, contradicting quasilocality of the Lindblad opera- 
tors. Thus, fik = ie^ with constant sign for all /c, corre- 
sponding to an average filling fik = |(1 — ^2,fc) = or 1. 
The situation is different for the second case: At the gap 
closing point, the vector fik points in the ±?/-direction 
for all /c, and there is no constraint on the purity. In 
fact, the system is half filled, fik = S = \i allowing 
thermodynamic observables to be properly defined. 
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